DirectionAngle Function

public function DirectionAngle(pointA, pointB) result(angle)

Uses

Find the direction angle (radians) between two points in a 2D space, measured in the clockwise direction. Formula gives values 0°<θ<90° for lines with positive slope and values 90°<θ<180° for lines with negative slope. The result does depend on which point is point 1 and which is point 2. Assume radians as input unit for geodetic reference system.

bearing angle on spheroid

Arguments

Type IntentOptional Attributes Name
type(Coordinate), intent(in) :: pointA

point A

type(Coordinate), intent(in) :: pointB

point B

Return Value real(kind=float)

bearing angle (radians)


Variables

Type Visibility Attributes Name Initial
real(kind=float), public :: a1
real(kind=float), public :: a2
real(kind=float), public :: b1
real(kind=float), public :: b2
real(kind=float), public :: deltaLongitude
real(kind=float), public :: x
real(kind=float), public :: y

Source Code

FUNCTION DirectionAngle &
!
(pointA, pointB) &
!
RESULT (angle)

USE Units, ONLY: &
!Imported parametes:
pi

IMPLICIT NONE

!Arguments with intent in:
TYPE (Coordinate), INTENT(IN) :: pointA !! point A 
TYPE (Coordinate), INTENT(IN) :: pointB !! point B 


!Local declarations:
REAL (KIND = float) :: angle !!bearing angle (radians)
REAL (KIND = float) :: a1, a2, b1, b2, x, y, deltaLongitude

!----------------------end of declarations-------------------------------------
a1 = pointA %  easting !longitude of point A
a2 = pointA %  northing !latitude of point A

b1 = pointB %  easting !longitude of point B
b2 = pointB %  northing !latitude of point B


IF ( pointA %  system % system == geodetic) THEN !!bearing angle on spheroid
    
    deltaLongitude = b1 - a1
    
    y = COS ( a2 ) * SIN (b2) - SIN ( a2 ) * COS ( b2 ) * COS ( deltaLongitude )
    
    x = SIN ( deltaLongitude ) * COS ( b2 )
    
    IF ( x == 0.) THEN
      angle = 0.
    ELSE
      angle = pi / 2. - ATAN ( y / x )
    END IF
    
ELSE  !points are in a projected coordinate system
    
    y = b2 - a2
    
    x = b1 - a1
    
    IF ( x == 0.) THEN
      angle = 0.
    ELSE
      angle = pi / 2. - ATAN ( y / x )
    END IF
    
END IF


RETURN
END FUNCTION DirectionAngle